#include "NeutrinoTools.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TDirectory.h"
#include "TSystem.h"

#include "style.h"

Bool_t ToReturn(const TString xvar, const TString cvar, Bool_t & kqacut0, Bool_t & kosci0)
{
  const TString pwd=gSystem->pwd();

  if(pwd.Contains("SK")){
    if(cvar=="PID"||cvar.Contains("Parent") || cvar.Contains("Nhits") || xvar.Contains("Nhits") || cvar.Contains("Flight") || xvar.Contains("Flight") || cvar.Contains("Charge") || xvar.Contains("Charge")  || cvar.Contains("LeptonType") || xvar.Contains("LeptonType")  ){
      return kTRUE;;
    }
  }

  if(pwd.Contains("SK")){
    kqacut0=kFALSE;
  }
  else{
    kosci0=kFALSE;
  }

  if(pwd.Contains("ND280")){
    if(!(cvar=="CCpurity" && (xvar=="fNeutMode" || xvar=="fNeutrinoType"))){
      if((cvar=="CCeff" || cvar=="CCpurity") && kqacut0==0){
        return kTRUE;;
      }
    }
  }

  if(cvar=="Q2PPt" && xvar=="fQ2Sim"){
    return kTRUE;;
  }

  if(cvar=="EnuPPt" && xvar=="fNeutrinoSim->E()"){
    return kTRUE;;
  }

//  if(cvar=="CCpurity" && (xvar=="fNeutMode" || xvar=="fNeutrinoType")){
//    return kTRUE;;
//  }

  if(cvar=="CCeff" && (xvar=="fNeutrinoType")){
    return kTRUE;;
  }

  if(cvar=="PID" && xvar=="f@TypeSim"){
    return kTRUE;;
  }

  return kFALSE;
}

TCanvas * getCanvas()
{
  TCanvas *c1=new TCanvas;
  style::SetGlobalStyle();
  style::PadSetup(c1);

  const Double_t lm=0.15;
  const Double_t rm=0.03;
  const Double_t tm=0.08;
  const Double_t bm=0.17;

  c1->SetLeftMargin(lm);
  c1->SetRightMargin(rm);
  c1->SetTopMargin(tm);
  c1->SetBottomMargin(bm);

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  gStyle->SetOptStat(0);

  c1->SetLogx();

  return c1;
}

void Config(Int_t & xnbin, Double_t & xmin, Double_t & xmax, Int_t & ncut, TString finalcut[], TString cname[], TString & header, const TString xvar0, const TString cvar, const Bool_t kqacut, const Bool_t kosci,const Bool_t kproton)
{
  TString xvar(xvar0);
  xvar.ReplaceAll("fMuon","f@");
  xvar.ReplaceAll("fProton","f@");

  xnbin = 50;
  xmin = 1e-3;
  
  const TString pwd=gSystem->pwd();
  Double_t ptmax = 3;
  if(pwd.Contains("ND280") && xvar.Contains("Rec")){
    ptmax = 3;//kqacut?3:20;
  }
  xmax = xvar.Contains("Pt")?ptmax:20;

  ncut = 0;
  TString cuts[10];
  TString cut0;
  header="";

  //const TString cutCCRec = Form("(f@NRec > 0 %s)", kqacut? Form(" && %s",(NeutrinoTools::GetQACut()).Data()): "");
  const TString cutCCRec = Form("(f@NRec > 0 %s)", kqacut? "badbad": "");
  const TString conCCRec = Form("@ rec %s",kqacut?"w/ QA":"w/o QA");

  //const TString cutCCpTRec = Form("(f@NRec > 0 %s)", kqacut? Form(" && %s",(NeutrinoTools::GetQACut()).Data()): Form(" && %s",(NeutrinoTools::GetPtResCut()).Data()));
  const TString cutCCpTRec = Form("(f@NRec > 0 %s)", "&& 1");
  const TString conCCpTRec = Form("@ rec %s",kqacut?"w/ p_{t}+PID QA":"w/ p_{t} QA");

  const TString cutCCSim = "( fNeutrinoType==14 && fNeutMode>0 && fNeutMode<30 && f@NSim>0 )";
  const TString conCCSim = "#nu_{#mu}CC";

  const TString cutIncl = "1";
  const TString conIncl = "";

  const TString cutCCFull = Form("%s && %s", cutCCSim.Data(), cutCCRec.Data());
  const TString conCCFull = conCCSim+conCCRec;

  const TString cutSimFull = xvar.Contains("Sim")? cutCCSim : cutCCFull;
  const TString conSimFull = xvar.Contains("Sim")? conCCSim : conCCFull;

  //================================================================================
  if(cvar == "CmpOsci" && pwd.Contains("SK")){
    ncut = 2;

    cuts[0]="fkOsciMuSurvive==1";
    cname[0]="#nu_{#mu}#rightarrow#nu_{#mu}";

    cuts[1]="fkOsciMuSurvive==0";
    cname[1]="#nu_{#mu}#rightarrow others";

    //everything
    cut0 = cutIncl;
    header=conIncl;
  }
  else if(cvar == "EnuNCCC"){
    ncut=2;

    cuts[0]= cutCCSim;
    cname[0]=conCCSim;

    cuts[1]= Form("!(%s)", cuts[0].Data());
    cname[1]=Form("#slash{%s}",cname[0].Data());

    xmin = 5e-2;

    //everything
    cut0 = cutIncl;
    header=conIncl;
  }
  //================================================================================
  else if(cvar == "CCeff"){
    ncut=2;

    cuts[0]= cutCCRec; 
    cname[0]=conCCRec;

    cuts[1]= Form("!(%s)", cuts[0].Data());
    cname[1]=Form("#slash{%s}",cname[0].Data());

    cut0 = cutCCSim;
    header=conCCSim;
  }
  //================================================================================
  else if(cvar == "CCpurity"){
    ncut=2;

    cuts[0]= cutCCSim;
    cname[0]=conCCSim;

    cuts[1]= Form("!(%s)", cuts[0].Data());
    cname[1]=Form("#slash{%s}",cname[0].Data());

    cut0 = cutCCRec;
    header = conCCRec;
  }
  //================================================================================
  else if(cvar == "PID" && pwd.Contains("ND280")){
    Int_t id = 0;

    cuts[id]="f@SimPDG==13";
    cname[id++]="#mu^{-}";
    
    cuts[id]="f@SimPDG==-13";
    cname[id++]="#mu^{+}";

    cuts[id]="abs(f@SimPDG)==11";
    cname[id++]="e";

    cuts[id]="abs(f@SimPDG)==211";
    cname[id++]="#pi";
    
    cuts[id]="abs(f@SimPDG)==2212";
    cname[id++]="p";
    
    ncut = id;

    //require also CCSim to avoid double counting with CCpurity
    cut0 = cutCCpTRec;
    header = conCCpTRec;
  }
  //================================================================================
  else if(cvar == "PtRes"){
    Int_t id = 0;
  
    cuts[id]="(f@RecPt->Mag()/f@SimPt->Mag())>=2";
    cname[id++]="p_{t}^{@ rec}/p_{t}^{@ sim}>2";
    
    cuts[id]="abs(f@RecPt->Mag()/f@SimPt->Mag()-1)>0.2 && (f@RecPt->Mag()/f@SimPt->Mag())<2";
    cname[id++]="<0.8, 1.2-2";
    
    cuts[id]="abs(f@RecPt->Mag()/f@SimPt->Mag()-1)<=0.2 && abs(f@RecPt->Mag()/f@SimPt->Mag()-1)>0.1";
    cname[id++]="0.8-0.9, 1.1-1.2";
    
    cuts[id]="abs(f@RecPt->Mag()/f@SimPt->Mag()-1)<=0.1";
    cname[id++]="0.9-1.1";

    ncut = id;

    cut0= cutCCFull;
    //header=conCCFull+" #frac{p_{t}^{@ rec}}{p_{t}^{@ sim}}";
    header=conCCFull;
  }
  else if(cvar == "PRes"){
    Int_t id = 0;
  
    cuts[id]="(f@Rec->P()/f@Sim->P())>=2";
    cname[id++]=">2";
    
    cuts[id]="abs(f@Rec->P()/f@Sim->P()-1)>0.2 && (f@Rec->P()/f@Sim->P())<2";
    cname[id++]="<0.8, 1.2-2";
    
    cuts[id]="abs(f@Rec->P()/f@Sim->P()-1)<=0.2 && abs(f@Rec->P()/f@Sim->P()-1)>0.1";
    cname[id++]="0.8-0.9, 1.1-1.2";
    
    cuts[id]="abs(f@Rec->P()/f@Sim->P()-1)<=0.1";
    cname[id++]="0.9-1.1";

    ncut = id;

    cut0= cutCCFull;
    header=conCCFull+" #frac{p_{@ rec}}{p_{@ sim}}";
  }
  else if(cvar == "MultiplicityFraction"){
    ncut=2;

    cuts[0]="fMultiplicity==1";
    cname[0]="N_{all rec}=1";

    cuts[1]="fMultiplicity>1";
    cname[1]="N_{all rec}>1";

    cut0 = cutCCFull;
    header=conCCFull;
  }
  else if(cvar == "MulPSelQ2" ){
    ncut=3;

    cuts[2]="fMultiplicity==1";
    cname[2]="N_{all rec}=1";

    cuts[1]="fMultiplicity>1 && f@Rec->P()<=1";
    cname[1]="N_{all rec}>1, p_{@ rec}#leq 1";

    cuts[0]="fMultiplicity>1 && f@Rec->P()>1";
    cname[0]="N_{all rec}>1, p_{@ rec}>1";

    cut0 = cutCCFull;
    header=conCCFull;
  }
  else if(cvar == "PtSelQ2"){
    ncut=2;

    cuts[1]="f@RecPt->Mag()<=1";
    cname[1]="p_{t}^{@ rec}#leq 1";

    cuts[0]="f@RecPt->Mag()>1";
    cname[0]="p_{t}^{@ rec}>1";

    cut0 = cutCCFull;
    header=conCCFull;
  }
  //================================================================================
  else if(cvar == "ModePPt"){
    Int_t id = 0;

    TString modename;

    cuts[id]=NeutrinoTools::GetModeCut(1,modename);
    cname[id++]=modename;

    cuts[id]=NeutrinoTools::GetModeCut(2,modename);
    cname[id++]=modename;

    cuts[id]=NeutrinoTools::GetModeCut(3,modename);
    cname[id++]=modename;

    cuts[id]=NeutrinoTools::GetModeCut(4,modename);
    cname[id++]=modename;

    cuts[id]=NeutrinoTools::GetModeCut(5,modename);
    cname[id++]=modename;

    //checked, no others!
    //cuts[id]=Form("!(%s) && !(%s) && !(%s) && !(%s)", cuts[0].Data(), cuts[1].Data(), cuts[2].Data(), cuts[3].Data());
    //cname[id++]="others";

    ncut=id;
    //for Rec, the NC component has mode >=30
    //for Sim, checked that NC ~ MuonNSim=0

    //cut0 =cutCCRec;
    //header = conCCRec;
    cut0 =cutSimFull;
    header = conSimFull;
  }
  else if(cvar == "ParentPPt" && pwd.Contains("ND280")){
    const Int_t parents[]={13, 130,211,321};
    const TString parname[]={"#mu", "K_{L}","#pi^{#pm}","K^{#pm}"};

    ncut = sizeof(parents)/sizeof(Int_t);

    for(Int_t icut=0; icut<ncut; icut++){
      cuts[icut]=Form("TMath::Abs(fNeutrinoParentPDG)==%d", parents[icut]);
      cname[icut]=parname[icut];
    }

    cut0 =cutSimFull;
    header = conSimFull;
  }
  else if(cvar == "Q2PPt" || cvar == "EnuPPt" ){
    TString tmpcv, tmpcn;
    if(cvar.Contains("Q2")){
      tmpcv = "Q^{2}";
      tmpcn = "fQ2Sim";
    }
    else if(cvar.Contains("Enu")){
      tmpcv = "E_{#nu}";
      tmpcn = "fNeutrinoSim->E()";
    }

    const Double_t nodes[]={-1e10, 0.2, 0.4, 1, 2, 1e10};
    ncut = sizeof(nodes)/sizeof(Double_t)-1;
    for(Int_t icut=0; icut<ncut; icut++){
      TString & cncn = cname[icut];
      TString & ctct = cuts[icut];
      if(fabs(nodes[icut])<10){
        cncn=Form("%.1f<%s", nodes[icut], tmpcv.Data());
      }
      else{
        cncn=Form("%s", tmpcv.Data());
      }

      if(fabs(nodes[icut+1])<10){
        cncn+=Form("<%.1f", nodes[icut+1]);
      }

      ctct=Form("%s>%f && %s<%f", tmpcn.Data(), nodes[icut], tmpcn.Data(), nodes[icut+1]);
    }


    cut0 =cutSimFull;
    header = conSimFull;
  }
  else{
    printf("wrong cvar %s\n", cvar.Data()); exit(1);
  }

  //add survival requirement, i.e. physics case at SK
  if(kosci && cvar!="CmpOsci"){
    //there is osci = -999 for ND
    cut0 += "&& fkOsciMuSurvive==1 ";
  }

  //===>
  if(xvar=="fNeutrinoSim->E()"){
    xmin = 0.99e-1;
  }
  
  if(xvar=="f@Rec->P()"){
    //xmin = 0.99e-1;
    xmin = 1e-3;
  }

  if(xvar=="fMultiplicity"){
    xmin = 0.5;
    
    if(pwd.Contains("SK")){
      xmax = 5.5;
    }
    else{
      xmax=10.5;
    }
    xnbin = xmax-xmin;
  }
  
  if(xvar=="fNeutMode"){
    xmin = -61.5;
    xmax =  61.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="fNeutrinoType"){
    xmin = -15.5;
    xmax =  15.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="(f@Charge[1]==-1)*2+(f@Charge[0]==-1)"){
    xmin = -0.5;
    xmax =  3.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="f@RecFlightPath->Pt()"){
    xmin = 0;
    xmax = 5;
    xnbin = 25;
  }

  if(xvar=="f@RecVertex->X()" || xvar=="f@RecVertex->Y()"){
    xmin = -1;
    xmax = 1;
    xnbin = 50;
  }

  if(xvar=="f@RecVertex->Z()"){
    xmin = 0.1;
    xmax = 0.5;
    xnbin = 50;
  }

  if(xvar=="f@RecEndPos->X()" || xvar=="f@RecEndPos->Y()"){
    xmin = -3;
    xmax = 3;
    xnbin = 60;
  }

  if(xvar=="f@RecEndPos->Z()"){
    xmin = -0.5;
    xmax = 4.5;
    xnbin = 80;
  }

  if(xvar=="f@RecFlightPath->Mag()"){
    xmin = 0;
    xmax = 6;
    xnbin = 30;
  }


  if(xvar=="fRangeMomentumMuon"){
    xmin = 0;
    xmax = 2;
    xnbin = 30;
  }

  if(xvar=="f@RecNhits[0]" || xvar=="fProtonRecNhits[0]"){
    xmin = 0;
    xmax = 300;
    xnbin=30;
  }

  if(xvar=="f@Length[0]"||xvar=="f@Length[1]"||xvar=="f@Length[2]"){
    xmin = -0.1;
    xmax = 3.1;
    xnbin=150;
  }

  if(xvar=="f@RecNhits[1]"){
    xmin = 0;
    xmax = 150;
    xnbin=75;
  }

  if(xvar=="f@TypeSim"){
    xmin = -6.5;
    xmax = 6.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="fNeutrinoParentType"){
    xmin = -6.5;
    xmax = 6.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="fNeutrinoRec->Theta()*TMath::RadToDeg()"){
    xmin = 1.3;
    xmax = 2.0;
    xnbin = 50;
  }

  if(xvar=="fNeutrinoSim->Theta()*TMath::RadToDeg()"){
    xmin = 0;
    xmax = 3;
    xnbin = 50;
  }

  if(xvar=="f@AlphaSim" || xvar=="f@AlphaRec" || xvar=="f@Rec->Theta()*TMath::RadToDeg()"){
    xmin = 0;
    xmax = 180;
    xnbin = 50;
  }

  if((cvar=="ModePPt" || cvar=="Q2PPt") ){
    if(xvar=="f@SimPt->Mag()"){
      xmin = 8e-2;
      xmax = 3;
      xnbin = 30;
    }
    else if(xvar=="f@Sim->P()"){
      xmin = 8e-2;
      xmax = 10;
      xnbin = 15;
    }
  }

  if(xvar.Contains("NDOF")){
    xmin = 1e-1;
    xmax = 10;
    xnbin = 20;
  }

  if(kproton){
    header.ReplaceAll("@","p");
  }
  else{
    header.ReplaceAll("@","#mu");
  }

  for(Int_t icut=0; icut<ncut; icut++){
    finalcut[icut]=(cut0+" && ")+cuts[icut];

    if(kproton){
      finalcut[icut].ReplaceAll("@","Proton");
      cname[icut].ReplaceAll("@","p");
    }
    else{
      finalcut[icut].ReplaceAll("@","Muon");
      cname[icut].ReplaceAll("@","#mu");
    }
  }
}

void plotStack(const TString fin, const TString xvar, const TString cvar, Bool_t kqacut0, const Bool_t kproton, Bool_t kosci0=kTRUE)
{
  if(ToReturn(xvar, cvar, kqacut0, kosci0)){
    return;
  }

  const Bool_t kqacut = kqacut0;
  const Bool_t kosci = kosci0;

  TFile::Open(fin);
  TTree * tree=(TTree*) gDirectory->Get("tree");

  TCanvas * c1 = getCanvas();

  //-----------------------------------------------------------
  Int_t xnbin;
  Double_t xmin;
  Double_t xmax;
  Int_t ncut;
  TString finalcut[10];
  TString cname[10];
  TString header;
  Config(xnbin, xmin, xmax, ncut, finalcut, cname, header, xvar, cvar, kqacut, kosci, kproton);
  
  //<===

  //check!
  //should only be a very small amount (~10) underflow
  //grep FLOW tmp | grep -v "UNDERFLOW 0 "  | grep UNDERFLOW
  //the non-0 is due to MuonN cut, and should all be the same value for the same cut: 3 possibilities: f@NSim>0, f@NSim>0 && f@NRec>0, f@NRec>0.
  //grep last tmp  | grep "d 0;" -v| grep " d "

  //----------------------------------------------------------------------
  const Int_t cols[]={kRed-4, kOrange, kCyan, kMagenta-7, kGreen};
  const Int_t mks[]={20,21,22,23,29};

  TLegend *lg = new TLegend(0.67, 0.72, 0.95, 0.90);
  style::ResetStyle(lg);
  lg->SetHeader(header);
  lg->SetBorderSize(0);
  lg->SetFillColor(0);

  TH1D *hcount[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0};

  //only for calculation of efficiency
  TH1D * htot = new TH1D("htot","", xnbin, xmin, xmax); 
  if(xvar!="fMultiplicity" && !xvar.Contains("RecNhits[0]") && !xvar.Contains("RecFlightPath->Pt()")){
    NeutrinoTools::BinLog(htot->GetXaxis());
  }
  style::ResetStyle(htot);

  const TString tottit(Form(";%s;Counts", (NeutrinoTools::GetTitleFromVar(xvar)).Data()));
  THStack * countStack = new THStack("countStack","");
  countStack->SetTitle(tottit);

  //some f@SimPt->Mag has ~ 1e-10 entries (~4), all with fNeutMode = 26, very parallel to neutrino  
  Double_t totn=0;
  for(Int_t ii=0; ii<ncut; ii++){
    hcount[ii] = new TH1D(Form("hh%d",ii),"", xnbin, xmin, xmax); 
    if(xvar!="fMultiplicity" && !xvar.Contains("RecNhits[0]") && !xvar.Contains("RecFlightPath->Pt()")){
      NeutrinoTools::BinLog(hcount[ii]->GetXaxis());
    }
    style::ResetStyle(hcount[ii]);

    printf("[%d]: xvar %s , finalcut { %s } , cname %s\n", ii, xvar.Data(), finalcut[ii].Data(), cname[ii].Data());
    tree->Draw(xvar+Form(">>%s",hcount[ii]->GetName()), finalcut[ii]);
   
    totn+=hcount[ii]->GetEntries();
    printf("\nxvar %s || cvar %s || finalcut { %s } || kqacut %d || %s %d: entries %.0f/%.0f d %.0f; UNDERFLOW %.0f OVERFLOW %.0f\n\n",xvar.Data(), cvar.Data(), finalcut[ii].Data(), kqacut,  ii==ncut-1? "last":"", ii, hcount[ii]->GetEntries(), (Double_t)(tree->GetEntries()), totn-(Double_t)(tree->GetEntries()), hcount[ii]->GetBinContent(0), hcount[ii]->GetBinContent(hcount[ii]->GetNbinsX()+1));

    hcount[ii]->SetFillColor(cols[(ii+1)%5]);
    hcount[ii]->SetLineColor(kBlack);
    //hcount[ii]->SetLineColor(cols[(ii+1)%5]);
    hcount[ii]->SetMarkerColor(cols[(ii+1)%5]); 
    hcount[ii]->SetMarkerStyle(mks[ii]);
    hcount[ii]->SetMarkerSize(2);
    hcount[ii]->SetLineWidth(1);

    hcount[ii]->SetTitle(cname[ii]+tottit);
    countStack->Add(hcount[ii]);
    lg->AddEntry(hcount[ii], cname[ii],"f");

    htot->Add(hcount[ii]);
  }

  //tmp
  countStack->Draw();
  //has to draw first
  style::ResetStyle(countStack);

  TString effTit(tottit);
  //effTit.ReplaceAll("Counts","#epsilon");
  effTit.ReplaceAll("Counts","Fraction");
  THStack * effStack = new THStack("effStack","");
  effStack->SetTitle(effTit);
  htot->Sumw2();
  for(Int_t ii=0; ii<ncut; ii++){
    printf("doing efficiency [%d]: %s , %s\n", ii, finalcut[ii].Data(), cname[ii].Data());
    TH1D * heff = (TH1D*) hcount[ii]->Clone(Form("%seff", hcount[ii]->GetName()));
    style::ResetStyle(heff);
    heff->Sumw2();
    heff->Divide(heff, htot, 1, 1, "B");
    heff->SetTitle(cname[ii]+effTit);
    effStack->Add(heff);
  }

  //-------------------------------------------------------------
  TString tmp(xvar);
  tmp.ReplaceAll("*TMath::RadToDeg()","");
  tmp.ReplaceAll(">","");
  tmp.ReplaceAll("(","");
  tmp.ReplaceAll(")","");
  tmp.ReplaceAll("-","");
  tmp.ReplaceAll("/","");
  tmp.ReplaceAll("[","");
  tmp.ReplaceAll("]","");
  tmp.ReplaceAll("+","");
  tmp.ReplaceAll("*","");
  tmp.ReplaceAll("=","");
  tmp.Prepend("outplot/");
  printf("tmp %s\n", tmp.Data());

  //c1->SetTopMargin(0.08);
 
  for(int ii=0; ii<2; ii++){
    gPad->SetLogy(ii); 
    countStack->SetMinimum(ii?0.9:0);
    for(Int_t jj=0; jj<2; jj++){
      gPad->SetLogx(jj);
      countStack->Draw("hist");
      lg->Draw();
      c1->Print(tmp+Form("countWithStack%sLogx%dLogy%dqacut%d%skproton%d.png",cvar.Data(), gPad->GetLogx(),  gPad->GetLogy(), kqacut,kosci?"Surv":"", kproton));
      
      countStack->Draw("nostack lpe");
      lg->Draw();
      c1->Print(tmp+Form("countNoStack%sLogx%dLogy%dqacut%d%skproton%d.png",cvar.Data(), gPad->GetLogx(),  gPad->GetLogy(), kqacut,kosci?"Surv":"", kproton));
    }
  }

  //tmp
  effStack->Draw();
  style::ResetStyle(effStack);

  //c1->SetTopMargin(0.02);

  for(int ii=0; ii<2; ii++){
    gPad->SetLogy(ii); 
    //effective, but the number is not accurate
    effStack->SetMinimum(ii?1e-3:0);
    effStack->SetMaximum(ii?15:1.41);
  
    for(Int_t jj=0; jj<2; jj++){
      gPad->SetLogx(jj);
      effStack->Draw("hist");
      lg->Draw();
      c1->Print(tmp+Form("effWithStack%sLogx%dLogy%dqacut%d%skproton%d.png",cvar.Data(), gPad->GetLogx(),  gPad->GetLogy(), kqacut,kosci?"Surv":"", kproton));
      
      effStack->Draw("nostack lpe");
      lg->Draw();
      c1->Print(tmp+Form("effNoStack%sLogx%dLogy%dqacut%d%skproton%d.png",cvar.Data(), gPad->GetLogx(),  gPad->GetLogy(), kqacut,kosci?"Surv":"", kproton));
    }
  }
}
